The aim of this tutorial is to show how to query and extract data, plug them in Pyhton tools and custom code
To install NDLab:
go to the GitHub page and read the How to install section
In some examples this Notebook makes use of two external packages: pandas and plotly.
You can of course use NDLab without them, they just make plotting easier.
To install plotly:
$ pip install plotly
$ pip install ipywidgets
$ jupyter labextension install jupyterlab-plotly
To install pandas:
$ pip install pandas
(You will have to restart jupyer lab)
Let's start importing NDLab:
# NuclearDataLab
import ndlab as nl
# this is needed only for using autocompletion when typing
from ndlaborm import *
import pandas as pd
import io
Just type the properties you want. As example, for the nuclide charge radius type NUCLIDE.CHARGE_RADIUS
Here the list of entites one can query
NUCLIDE, LEVEL, GAMMA, L_DECAY, CUM_FY, IND_FY
plus the set of decay radiations:
DR_ALPHA, DR_GAMMA, DR_BETAM, DR_BETAP, DR_ANTI_NU, DR_NU, DR_X, DR_ANNHIL, DR_AUGER, DR_CONV_EL, DR_DELAYED, DR_PHOTON_TOTAL
In the picture below the fields providing the links are in green - to nuclide, and red - to level.

print(description(NUCLIDE))
NUCLIDE: properties of the nuclide in its ground state ABUNDANCE *(q)* natural abundance, in mole fraction ATOMIC_MASS *(Q)* atomic mass, in micro AMU BETA_DECAY_EN *(Q)* energy available for beta decay [keV] BINDING_EN *(Q)* binding energy per nucleon [keV] CHARGE_RADIUS *(Q)* Root-mean-square of the nuclear charge radius, expressed in fm. ELEM_SYMBOL *(S)* symbol of the element MASS_EXCESS *(Q)* mass excess [keV] N number of neutrons NUC_ID *(S)* identifier, e.g. 135XE QA *(Q)* alpha decay Q energy [keV] QBMN *(Q)* beta- + n decay energy [keV] QEC *(Q)* electron capture Q-value [keV] S2N *(Q)* 2-neutron separation energy [keV] S2P *(Q)* 2-proton separation energy [keV] SN *(Q)* neutron separation energy [keV] SP *(Q)* proton separation energy [keV] Z number of protons
The simplest way to retrieve data is to call the nl.csv_data(fields, filter) function, or the nl.json_data(fields, filter) function
<div>
fields = "GAMMA.NUC.Z as z, GAMMA.NUC.N as n , abs( GAMMA.MIXING_RATIO / GAMMA.ENERGY / 0.835 * 1000) as m"
filter = "( GAMMA.NUC.Z % 2 = 0 ) and ( GAMMA.NUC.N % 2 = 0 ) " # gamma from even-even nuclides
filter += " and GAMMA.START_LEVEL.JP = '2+' and GAMMA.START_LEVEL.JP_ORDER = 2 "# starts from Jp = 2+ , 2nd occurrence
#filter += # ends at Jp = 2+ , 1st occurrence
#filter += # E2 + M1 multipolarity
#print(nl.csv_data(fields, filter)[0:250])
So far only
the NDLab package has been used, this is enough to retrieve data and export them in CSV or Json
What follows shows how to embed NDLab in Pandas and Plotly
Run the cell below which loads the packages
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
# Plotly
import plotly.graph_objects as go
import plotly.express as px
The Y axis of the Krane's picture above needs Log(Δ).
Using pandas, numpy, and matplotlib data are prepared and plotted.
With fields and filter of the above practice, the following cell loads the data and shows part of the dataframe
fields = " GAMMA.NUC.Z as z, GAMMA.NUC.N as n, abs( GAMMA.MIXING_RATIO / GAMMA.ENERGY / 0.835 * 1000 ) as y " # z, n, mixing / energy
filter = "(GAMMA.NUC.Z % 2 = 0 ) AND (GAMMA.NUC.N % 2 = 0 ) " # even-even nuclides
filter += " and GAMMA.START_LEVEL.JP = '2+' and GAMMA.START_LEVEL.JP_ORDER = 2 "# starts from Jp = 2+ , 2nd occurrence
filter += " and GAMMA.END_LEVEL.JP = '2+' and GAMMA.END_LEVEL.JP_ORDER = 1 " # ends at Jp = 2+ , 1st occurrence
filter += " and ( GAMMA.MULTIPOLARITY = 'E2+M1' or GAMMA.MULTIPOLARITY = 'M1+E2' )" # E2 + M1 multipolarity
# load the data into a dataframe
df = nl.pandas_df(fields, filter, pd)
df
| z | n | y | |
|---|---|---|---|
| 0 | 42 | 58 | 9.975355 |
| 1 | 44 | 56 | 5.386445 |
| 2 | 46 | 56 | 3.429497 |
| 3 | 46 | 64 | 12.527247 |
| 4 | 48 | 64 | 6.893959 |
| ... | ... | ... | ... |
| 129 | 40 | 50 | 0.266848 |
| 130 | 42 | 52 | 2.411730 |
| 131 | 42 | 54 | 0.732317 |
| 132 | 38 | 58 | 3.461286 |
| 133 | 40 | 56 | 0.226757 |
134 rows × 3 columns
Now the data are plotted by calculating Z + N for the x, and log(Δ) the the y
plt.scatter( df["z"] + df["n"], np.log(df["y"])*10, color='red' )
<matplotlib.collections.PathCollection at 0x7fe7d910f730>
It is worth to superimpose the new data with the old ones
plt.subplots(figsize=(28, 10))
plt.scatter( df["z"] + df["n"], np.log(df["y"])*10, color='red')
plt.imshow(mpimg.imread("../docs/_static/krane_2.png"), extent=[38, 205, -63, 79])
<matplotlib.image.AxesImage at 0x7fe8016a4700>
The paper claims that there are minima at closed shells
On a 2D plot it is difficult to see where the closed shells in N or Z are placed
The cell below uses plotly to plot a 3D picture where nuclides with closed shells are in red:
Plotly graphs are interactive, can be rotated and zoomed
# condition for closed-shell
filter2 = filter + " and ( GAMMA.NUC.Z in ( 2,8,20,28,50,82,126 ) or GAMMA.NUC.N in ( 2,8,20,28,50,82,126 ) )"
# load data
df2 = nl.pandas_df(fields, filter2, pd)
# plot 3D
fig = go.Figure(data=[
# all nuclides
go.Scatter3d(
x=df["z"], y=df["n"], z=np.log10(df.y)*10,
mode='markers', marker=dict(size=8, color='blue')
)
,
# the magic ones in red
go.Scatter3d(
x=df2["z"], y=df2["n"], z=np.log10(df2.y)*10,
mode='markers',marker=dict(size=8,color="red")
)
])
# set layout properties
fig.update_layout(margin=dict(l=0, r=0, b=30, t=0) ,width=800, height=800,
scene = dict(xaxis_title='Z',yaxis_title='N',zaxis_title='Delta')
).show()
Once you have a dataframe, you can export into almost anything
df.to_csv()
df.to_json()
df.to_hdf()
df.to_excel()
df.to_sql()
df.to_xml()
fields = "GAMMA.ENERGY as energy"
condition = "GAMMA.NUC.Z < 80 and GAMMA.NUC.Z > 0 and GAMMA.ENERGY < 2000"
dfg = nl.pandas_df(fields,condition,pd)
px.histogram( dfg, x="energy", nbins=200).show()
The picture above relies on your brain to find patterns and outliers, whilst
the cell below shows how to automatise the process.
The assumption is that the curve should be smooth, meaning that the first dervative should not have sudden jumps.
Comments in cell explain the steps
fields = "GAMMA.ENERGY as energy"
condition = "GAMMA.ENERGY BETWEEN 0 AND 2000"
# load a pandas dataframe using the convenience funtion defined above
a = nl.pandas_df(fields,condition,pd)
# bin the data (consult the web for further examples and explanations)
dft = pd.cut(a[a.columns[0]], bins=300).apply(lambda x: x.mid).value_counts(sort=False).to_frame()
dft = dft.rename(columns={"energy": "counts"})
# get the derivative
dft["delta"] = abs(dft.diff(periods = 1).counts )
# visualise the outlier (consult the web to check Plotly tools for data processing)
fig = px.box(dft, y=dft.delta).show()
# print the energy where the outlier occurs
print("Outlier energy ", dft.index[dft.index.get_loc(dft.delta.idxmax()) - 1 ])
Outlier energy 510.00550000000004
In the plot the transition type is not decoded: 1NU stands for 1st non-unique, etc...
fields = "DR_BETA.TRANS_TYPE as t , DR_BETA.LOGFT as l"
filter = " DR_BETA.LOGFT != 0" # this is just to avoid empty values
df = nl.pandas_df(fields, filter,pd)
# group by transition type
ans = [pd.DataFrame(y) for x, y in df.groupby(df.columns[0], as_index=False)]
# plot log_ft for each transition type
fig = go.Figure()
for a in ans:
tst = pd.cut(a['l'], bins=120).apply(lambda x: x.mid).value_counts(sort=False)
fig.add_trace(go.Scatter(x=tst.index.values, y = tst, mode="markers",name=a.iloc[0, 0]))
fig.update_layout(width=800, height=400, xaxis_range=[3,15] ).show()
fields = " DR_BETA.LOGFT as log_ft"
filter = "( DR_BETA.PARENT_LEVEL.JP = '0+') and DR_BETA.TRANS_TYPE ='A'"
df = nl.pandas_df(fields, filter, pd)
px.histogram(df, x=df.log_ft, nbins=80).update_layout(width=350, height=350).show()
With pandas it is easy to dump the data in almost any format. In this way one is then free to use them as input for other tools
Check the web for how to use the various writing - reading methods listed below
filter = ("LEVEL.NUC.N % 2 = 0 AND LEVEL.NUC.Z % 2 = 0 AND LEVEL.SEQNO = 1 AND LEVEL.JP != '2+' AND LEVEL.JP_METHOD = JP_STRONG ")
df = nl.pandas_df_nl(nl.levels(filter), pd)
df
| z | n | nucid | l_seqno | energy | energy_unc | energy_limit | half_life | half_life_unc | half_life_limit | ... | jp_str | quadrupole_em | quadrupole_em_unc | quadrupole_em_limit | dipole_mm | dipole_mm_unc | dipole_mm_limit | questionable | configuration | isospin | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 36 | 64 | 100KR | 1 | 309.00 | 10.0000 | = | 0.00 | 0.00 | = | ... | (2+) | 0.0 | 0.0 | = | 0.0 | 0.0 | = | NaN | NaN | NaN |
| 1 | 38 | 62 | 100SR | 1 | 129.18 | 0.0009 | = | 3.91 | 0.16 | = | ... | (2+) | 0.0 | 0.0 | = | 0.0 | 0.0 | = | NaN | NaN | NaN |
| 2 | 38 | 64 | 102SR | 1 | 126.00 | 0.2000 | = | 3.00 | 1.20 | = | ... | (2+) | 0.0 | 0.0 | = | 0.0 | 0.0 | = | NaN | NaN | NaN |
| 3 | 40 | 64 | 104ZR | 1 | 139.30 | 0.0300 | = | 2.00 | 0.30 | = | ... | (2+) | 0.0 | 0.0 | = | 0.0 | 0.0 | = | NaN | NaN | NaN |
| 4 | 52 | 54 | 106TE | 1 | 664.80 | 0.0300 | = | 0.00 | 0.00 | = | ... | (2+) | 0.0 | 0.0 | = | 0.0 | 0.0 | = | Y | NaN | NaN |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 155 | 74 | 118 | 192W | 1 | 219.00 | 0.0000 | = | 0.00 | 0.00 | = | ... | [2+] | 0.0 | 0.0 | = | 0.0 | 0.0 | = | NaN | NaN | NaN |
| 156 | 60 | 100 | 160ND | 1 | 65.20 | 0.0500 | = | 0.00 | 0.00 | = | ... | (2+) | 0.0 | 0.0 | = | 0.0 | 0.0 | = | NaN | NaN | NaN |
| 157 | 62 | 102 | 164SM | 1 | 69.00 | 1.0000 | = | 0.00 | 0.00 | = | ... | (2+) | 0.0 | 0.0 | = | 0.0 | 0.0 | = | NaN | NaN | NaN |
| 158 | 92 | 124 | 216U | 1 | 2240.00 | 42.0000 | = | 0.70 | 0.80 | = | ... | (8+) | 0.0 | 0.0 | = | 0.0 | 0.0 | = | NaN | NaN | NaN |
| 159 | 18 | 12 | 30AR | 1 | 700.00 | 0.0000 | = | 0.00 | 0.00 | = | ... | (2+) | 0.0 | 0.0 | = | 0.0 | 0.0 | = | NaN | NaN | NaN |
160 rows × 28 columns
Often one wants to write had hoc processing code that differs from filtering, grouping, binning, or any other manipulation. For example, calculate the energy balance of a decay.
Ndlab provides as a set of python classes and functions pushing the data extraction to the background. The functions accept the filter parameter (default no filter) to restrict the selection.
delayed_n = nl.dr_delayeds( "DR_DELAYED.TYPE = DELAY_N ")
parent_nucs = [ ld.parent for ld in delayed_n]
parent_nucs = nl.remove_doublers(parent_nucs)
The following is a list of hands-on examples. To see the full list of classes and function, please consult the guide
With NDLab's set of classes and functions, you do not need to take care of the data retrieval. See how you can calculate the energy balance of a decay.
filter = ""
# load a nuclide using its identifier
nuc = nl.nuclide("82SR")
# the decays() function loads an array with all the decays. Take the first one
decay = nuc.decays[0]
# all possible radiation types for the chosen decay
rads = [
decay.xs(), # X-ray
decay.gammas(), # Gamma
decay.convels(), # Conversion electron
decay.augers(), # Auger
decay.alphas(), # Alpha
decay.betas_m(), # B-
decay.betas_p(), # B+
decay.anti_nus(), # anti-neutrino
decay.nus(), # neutrino
decay.annihil() # annihilation
]
# measured energy for each radiation type
# note that for neutrini there are two contributions:
# - one for the beta+ process, included in the first term below
# - one for the electron capture process, included in the second
rads_en = [(sum([r[i].energy * r[i].intensity for i in range(len(r))])) for r in rads] + [(nu.energy_ec * nu.intensity_ec) for nu in decay.nus()]
# total measured energy
tot_en = sum(rads_en) / 100
# intensities are per 100 decay of the parent, need to renormalise Q with branching ratio
q_br = decay.q_togs * decay.perc/100
# difference with deposited energy
delta = ((q_br - tot_en)/q_br)*100
print("Branching Ratio: ", decay.perc/100 )
print("Energy accounted for [keV]" , tot_en )
print("Qb- * B.R. [keV] " , q_br)
print("Delta [%] " , int(delta.value * 100) /100)
print("\nNotice the built-in uncertainty propagation")
Branching Ratio: 1.0+/-0 Energy accounted for [keV] 13.31+/-0.05 Qb- * B.R. [keV] 178+/-7 Delta [%] 92.52 Notice the built-in uncertainty propagation
Summary of decay chain, see the next cell for a more sophisticated example
# the daughters and parents of a nuclide are lists with the direct daughters and parents, respectively
# the list contains the Ground States of the daughters/parents
nuc = nl.nuclide("238U")
print('Daughters and H-l [s]')
[print(d.nucid, d.gs.half_life_sec) for d in nuc.daughters]
print('\nParents and H-l [s]')
[print(p.nucid, p.gs.half_life_sec) for p in nuc.parents]
# the daughters_chain and parents_chain of a nuclide are lists with all daughters and parents, respectively
# the list contains the Ground States of the daughters/parents
print('\nDaughters chain \nNote that U-238 is listed as a daughter since it has an isomeric state that decays to the GS ')
dau = [n.gs for n in nl.nuclide("238U").daughters_chain]
for d in dau:
print(d.nucid)
print('\nParents chain')
par = [n.gs for n in nl.nuclide("205TL").parents_chain]
for p in par:
print(p.nucid)
Daughters and H-l [s] 234TH (2.0822+/-0.0026)e+06 238U (1.4100+/-0.0019)e+17 Parents and H-l [s] 238PA 137+/-6 238U (1.4100+/-0.0019)e+17 242PU (1.183+/-0.006)e+13 Daughters chain Note that U-238 is listed as a daughter since it has an isomeric state that decays to the GS 234TH 238U 234PA 234U 230TH 226RA 222RN 218PO 214PB 218AT 214BI 210TL 214PO 210PB 209PB 206HG 210BI 206TL 206PB 210PO 209BI 205TL 218RN Parents chain 205HG 205PB 209BI 205AU 206PT 205PT 205IR 205BI 209PO 205PO 209AT 205AT 209RN 205RN 209FR 209RA 205FR 213TH 217U 209AC 213PA 213AC 217PA 213RA 217TH 221U 213FR 217AC 221PA 225NP 229AM 233BK 213RN 217RA 221TH 225U 229PU 233CM 237CF 241FM 209PB 213AT 223RA 210TL 209TL 213PO 223AC 223FR 227TH 227PA 231NP 231PU 235AM 239BK 239CF 243ES 243FM 247MD 223RN 227AC 223AT 223PO 223BI 224BI 227RA 231PA 231TH 227FR 231AC 235U 231RA 231FR 235PA 235NP 239PU 235TH 235PU 239AM 239CM 243CF 247FM 251NO 255RF 259SG 263HS 267DS 243BK 247ES 251MD 255LR 259DB 239NP 243CM 239U 243AM 239PA 243PU 247BK 244AM 243NP 247CM 243U 247AM 251CF 247PU 251BK 251ES 255FM 251CM 255ES 259MD 255CF 259NO 263RF 251FM 255MD 255NO 259RF 263SG 267HS 271DS 259LR 263DB 267BH 247CF 227RN 227AT 227PO 231U 211HG 214BI 214PB 218AT 215TL 214TL 218PO 216HG 215HG 214HG 218BI 222RN 218PB 222AT 226RA 226FR 226AC 230TH 226RN 226AT 226PO 230PA 230AC 234U 230RA 230FR 234PA 234NP 238PU 234TH 234AC 238U 234RA 238PA 242PU 238TH 242NP 242AM 246CM 242U 246AM 246BK 250CF 246PU 250CM 254CF 254ES 258MD 262LR 266DB 270BH 274MT 278RG 282NH 246CF 250ES 246ES 250FM 246FM 250MD 246MD 250NO 254RF 258SG 254LR 258DB 262BH 266MT 254NO 258RF 262SG 266HS 270DS 250BK 254FM 254MD 258LR 262DB 266BH 270MT 274RG 278NH 234PU 238AM 234AM 238CM 234CM 238CF 234BK 242CF 238BK 242ES 238NP 242CM 242BK 213BI 209HG 217AT 213PB 217PO 221FR 217BI 221RN 217PB 217TL 221AT 225RA 221PO 221BI 222BI 225FR 229TH 225RN 225AT 229RA 225PO 229FR 229RN 229AT 229AC 229PA 233U 229U 233NP 229NP 233PU 233AM 237CM 241CF 241ES 245FM 245MD 237AM 241BK 245ES 249MD 253LR 257DB 261BH 233PA 237PU 237NP 233TH 237U 241AM 237PA 241PU 237TH 241NP 245CM 241U 245AM 245BK 249CF 245PU 249BK 249CM 253ES 253CF 257FM 257ES 257MD 257NO 257LR 261RF 261DB 257RF 261SG 265HS 269DS 265SG 269HS 273DS 277CN 253FM 253MD 253NO 245CF 249ES 249FM 241CM 233AC 233RA 233FR 225AC 225TH 213TL 213HG 210AU 209AU 217RN 221RA 217FR 221AC 225PA
The function decay_chain(nucid, level, initial population), at the bottom of the next cell, produces a description of a decay chain.
The relevant entities are saved in a data structure for further analysis, see the final printout.
Note that in the Th-234 step, one needs to consider the Pa-234 metastable feeding
Indeed the Th-234 -> Pa-234 decay has two lines (#2 and #15), as well as the Pa-234 -> U-234 (#3 and #16)
This example shows again how NDLab allows to focus on the Physics of the issue, and not on libraries' formats, structure, conventions ...
# dictionary filled by the decay_chain() function to keep the relevant entities for further analysis
SAVE = []
# Keeps track of the decays already processed
DONE = []
# the H-l threshold for considering a metastable, in seconds
H_L_META = "60"
# do not consider dacays below thisthershold, in %
DECAY_PERC_TRESHOLD = 5
# filter any alpha, beta the decay radiation feeding a metastable : the seqno of the fed level must not be 0, and its H-l must be above the threshold
FILTER_META = "( DECAY_RAD.DAUGHTER_FED_LEVEL.HALF_LIFE_SEC > " + H_L_META + " and DECAY_RAD.DAUGHTER_FED_LEVEL_SEQNO != 0 )"
# same as above but for gamma from IT decays
FILTER_META_GAMMA = "( DR_GAMMA.END_LEVEL.HALF_LIFE_SEC > " + H_L_META + " and DR_GAMMA.END_LEVEL_SEQNO != 0 )"
# keeps track of nuclides already populated, to renormalise the production rate
ALREADY_EXISTING = {}
# whether to print the decay balance
PRINT_BALANCE = True
# whether any radiation populates a metastable, consider only alpha, beta +/-, and isomeric transition
def get_rads_to_meta(d):
if(d.code in [DECAY_A]):
rads = d.alphas(FILTER_META)
elif (d.code in [DECAY_Bm]):
rads = d.betas_m(FILTER_META)
elif (d.code in [DECAY_Bp]):
rads = d.betas_p(FILTER_META)
elif (d.code in [DECAY_IT]):
rads = d.gammas(FILTER_META_GAMMA)
else:
print(' wrong decay ',d.code) # decay modes like delayed emission are not considered
rads = []
if(len(rads) > 1): print(rads[0].parent_nucid) # allow only 1 meta , warning if more
return rads
# the gamma feeding of a level from the above levels. For metastable, this intensity needs to be added to the direct feeding
def gamma_feed(decay, end_level):
filter_gm = "DR_GAMMA.END_LEVEL_SEQNO = " + str( end_level)
feeding = sum( g.intensity for g in decay.gammas(filter_gm))
return feeding
# energy balance of the decay, to see whether the decay schema is complete
def energy_balance(decay):
tot_en = decay.tot_measured_en()
if(tot_en == 0):
return False
print('\n', 'Decay', decay.mode.name,'from', decay.nucid, 'level ',decay.l_seqno, 'to', decay.daughter_nucid , str(decay.perc)+' %')#, ' cumulative %',str((cum_intensity*100))+' %' )
if( not PRINT_BALANCE): return True
# intensities are per 100 decay of the parent, need to renormalise Q with branching ratio
q_br = decay.q_togs * decay.perc/100
# difference with deposited energy
delta = ((q_br - tot_en)/q_br)*100
print("Branching Ratio: ", decay.perc/100 )
print("Energy accounted for [keV]" , tot_en )
print("Qb- * B.R. [keV] " , q_br)
print("Delta [%] " , int(delta.value * 100) /100)
return True
# this performs the business
# it is called recursively from parent to daughter(s). A daughter can be processed more than once if there are many pathways to it
# nucid , l_seqno are from the daughter of the previous decay, that now becomes the parent
# cum_intensity is the cumulative intensity, normalised to 100, with which the nuclide is produced:
def decay_chain(nucid, l_seqno, cum_intensity):
# check if the cumulative intensity has a contribution from another decay already processed
if nucid in ALREADY_EXISTING:
cum_intensity = ALREADY_EXISTING[nucid]
cum_intensity = (cum_intensity/100)
# get all decays
nuc = nl.nuclide(nucid)
decays = nuc.decays
if( len(decays)==0) : return
decays_todo = []
# process only the decays from the given level
for d in decays:
if(d.perc > DECAY_PERC_TRESHOLD and d.l_seqno == l_seqno): decays_todo.append(d)
for d in decays_todo:
# skip if the decay was already processed
if(d.pk in DONE or not energy_balance(d)): return
print('Cumulative population of the father:',str((cum_intensity*100))+' %' )
# get the decay to a metastable
rads_meta = get_rads_to_meta(d)
# total intesitiy that populates the metastable : direct feeding plus gamma from above
tot_meta_feed = 0
for r in rads_meta:
# no metastable populated, skip
if(len(rads_meta) == 0): break
tot_meta_feed = gamma_feed(d,r.fed_level.l_seqno) + r.intensity
print(' metastable fed from parent level', r.parent_l_seqno, ' --> to daughter level ',r.fed_level.l_seqno, ' gamma feeding: ' ,gamma_feed(d,r.fed_level.l_seqno) ,' direct feeding: ', r.intensity)
# check if meta and the GS have the same daughter
# in case, save the total BR for the daughter: when the daughter is processed it has the correct cumulative intensity
dau_decays = nl.nuclide(d.daughter_nucid).decays
dau_dau = ''
for dd in dau_decays:
if(dd.l_seqno == 0): # decay from GS of the daughter
dau_dau = dd.daughter_nucid
if(dd.l_seqno == r.fed_level.l_seqno and dd.daughter_nucid == dau_dau): # daughter of meta = daughter of GS
# simplify assuming that the meta has 2 decay modes, one IT to the GS, and one to the the same daughter of the GS
# example 234PA. More complex cases are not analysed here
ALREADY_EXISTING[dau_dau] = d.perc * cum_intensity
# remember that the decay has been processed
DONE.append(d.pk)
SAVE.append({'parent': nuc, 'level': l_seqno, 'fed_level':r.fed_level.l_seqno,'cumulative':(tot_meta_feed * cum_intensity), 'decay': d})
# go to the next step in the chain using the decay from the meta
decay_chain(d.daughter_nucid, r.fed_level.l_seqno, tot_meta_feed * cum_intensity)
# for the GS take the decay BR minus what goes to a meta
tot_gs_feed = d.perc - tot_meta_feed
# remember that the decay has been processed
DONE.append(d.pk)
SAVE.append({'parent': nuc, 'level': l_seqno, 'fed_level':0,'cumulative':(tot_gs_feed * cum_intensity), 'decay': d})
# go to the next step in the chain using the decay from GS
decay_chain(d.daughter_nucid, 0,tot_gs_feed * cum_intensity)
###########################
# CALL THE FUNCTION
# 238U is the initial nuclide,
# 0 is its decaying level seqno (0 = Ground State)
# 100 is the initial normalisation of the number of nuclides
decay_chain('238U',0, 100)
print('\nThe Dictionary has saved the relevant entities:\n')
cnt = 1
for s in SAVE:
print('*'+str(cnt),s['parent'].nucid, s['decay'].daughter_nucid, s['level'], s['fed_level'], s['cumulative'])
cnt = cnt+1
Decay A from 238U level 0 to 234TH 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] (4.20+/-0.18)e+03
Qb- * B.R. [keV] 4269.9+/-2.1
Delta [%] 1.64
Cumulative population of the father: 100.0 %
Decay B- from 234TH level 0 to 234PA 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] 66.8+/-1.6
Qb- * B.R. [keV] 274.0+/-3.0
Delta [%] 75.63
Cumulative population of the father: 100.0+/-0 %
metastable fed from parent level 0 --> to daughter level 2 gamma feeding: 4.23+/-0.28 direct feeding: 78.0+/-2.0
Decay B- from 234PA level 2 to 234U 99.84+/-0.04 %
Branching Ratio: 0.9984+/-0.0004
Energy accounted for [keV] ? 1656.4+/-3.2
Qb- * B.R. [keV] 2264+/-4
Delta [%] 26.84
Cumulative population of the father: 82.2+/-2.0 %
Decay A from 234U level 0 to 230TH 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] ? 4773+/-9
Qb- * B.R. [keV] 4857.5+/-0.7
Delta [%] 1.73
Cumulative population of the father: 100.0+/-0 %
Decay A from 230TH level 0 to 226RA 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] ? 4679+/-19
Qb- * B.R. [keV] 4770.0+/-1.5
Delta [%] 1.9
Cumulative population of the father: 100.0+/-0 %
Decay A from 226RA level 0 to 222RN 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] 4785+/-5
Qb- * B.R. [keV] 4870.70+/-0.25
Delta [%] 1.76
Cumulative population of the father: 100.0+/-0 %
Decay A from 222RN level 0 to 218PO 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] ? 5485.1+/-0.5
Qb- * B.R. [keV] 5590.40+/-0.30
Delta [%] 1.88
Cumulative population of the father: 100.0+/-0 %
Decay A from 218PO level 0 to 214PB 99.980+/-0.020 %
Branching Ratio: 0.99980+/-0.00020
Energy accounted for [keV] 6001.34+/-0.16
Qb- * B.R. [keV] 6113.5+/-1.2
Delta [%] 1.83
Cumulative population of the father: 100.0+/-0 %
Decay B- from 214PB level 0 to 214BI 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] 1016+/-10
Qb- * B.R. [keV] 1018+/-11
Delta [%] 0.15
Cumulative population of the father: 99.980+/-0.020 %
Decay B- from 214BI level 0 to 214PO 99.979+/-0.013 %
Branching Ratio: 0.99979+/-0.00013
Energy accounted for [keV] ? 3255+/-11
Qb- * B.R. [keV] 3268+/-11
Delta [%] 0.41
Cumulative population of the father: 99.980+/-0.020 %
Decay A from 214PO level 0 to 210PB 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] 7686.82+/-0.09
Qb- * B.R. [keV] 7833.54+/-0.06
Delta [%] 1.87
Cumulative population of the father: 99.959+/-0.024 %
Decay B- from 210PB level 0 to 210BI 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] 58.1+/-1.6
Qb- * B.R. [keV] 63.5+/-0.5
Delta [%] 8.54
Cumulative population of the father: 99.959+/-0.024 %
Decay B- from 210BI level 0 to 210PO 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] 1161+/-4
Qb- * B.R. [keV] 1161.2+/-0.8
Delta [%] -0.01
Cumulative population of the father: 99.959+/-0.024 %
Decay A from 210PO level 0 to 206PB 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] 5304.39+/-0.07
Qb- * B.R. [keV] 5407.53+/-0.07
Delta [%] 1.9
Cumulative population of the father: 99.959+/-0.024 %
Decay B- from 234PA level 0 to 234U 100.0+/-0 %
Branching Ratio: 1.0+/-0
Energy accounted for [keV] ? (2.53+/-0.04)e+03
Qb- * B.R. [keV] 2194+/-4
Delta [%] -15.48
Cumulative population of the father: 17.8+/-2.0 %
The Dictionary has saved the relevant entities:
*1 238U 234TH 0 0 100.0+/-0
*2 234TH 234PA 0 2 82.2+/-2.0
*3 234PA 234U 2 0 82.1+/-2.0
*4 234U 230TH 0 0 100.0+/-0
*5 230TH 226RA 0 0 100.0+/-0
*6 226RA 222RN 0 0 100.0+/-0
*7 222RN 218PO 0 0 100.0+/-0
*8 218PO 214PB 0 0 99.980+/-0.020
*9 214PB 214BI 0 0 99.980+/-0.020
*10 214BI 214PO 0 0 99.959+/-0.024
*11 214PO 210PB 0 0 99.959+/-0.024
*12 210PB 210BI 0 0 99.959+/-0.024
*13 210BI 210PO 0 0 99.959+/-0.024
*14 210PO 206PB 0 0 99.959+/-0.024
*15 234TH 234PA 0 0 17.8+/-2.0
*16 234PA 234U 0 0 17.8+/-2.0
Let's say that the level #15 of Xe-135 is populated by some reaction.
The following code gives the list of the levels populated by the gamma cascade, with energy and intensity
def cascade( my_level):
global todo
global levels
global done
for g in my_level.gammas() :
if(g.end_level.l_seqno != 0 and not (g.end_level.l_seqno in todo)):
todo[g.end_level.l_seqno] = g.end_level
if not g.end_level.l_seqno in result:
result[g.end_level.l_seqno] = result[my_level.l_seqno] * g.rel_photon_intens/100
else:
result[g.end_level.l_seqno] += result[my_level.l_seqno] * g.rel_photon_intens/100
if not my_level.l_seqno in done:
done.append(my_level.l_seqno)
todo = dict(sorted(todo.items(), reverse=True))
myrun = todo.copy()
for r in myrun:
if( not r in done):
cascade(levels[r])
return
nuc_id = "135XE"
start_level = 15
levels = nl.nuclide(nuc_id).levels()
result = {start_level : 1.0}
todo = {start_level : levels[start_level] }
done = []
cascade(levels[start_level])
# just printing
line = '|'.join(str(x).ljust(24) for x in [ "energy [keV] ", " population %", "level #"])
print(line +'\n')
for k in sorted(result):
line = '|'.join(str(x).ljust(24) for x in [str(levels[k].energy ), str(result[k]*100.0),levels[k].pk])
print(line)
energy [keV] | population % |level # 0.0+/-0 | 100+/-8 |135XE-0 288.455000+/-0.000015 | 0.0041+/-0.0007 |135XE-1 526.551000+/-0.000013 | 34+/-5 |135XE-2 1131.512000+/-0.000011 | 7+/-5 |135XE-3 1260.416000+/-0.000013 | 0.137+/-0.024 |135XE-4 1565.288000+/-0.000016 | 37+/-5 |135XE-8 1927.294000+/-0.000023 |100.0 |135XE-15
Nuclides often have many decay modes, either because of metastable states or because of a single energy state has more branchings
In gamma spectroscopy it is useful to have the total intensity of an energy line in the decay of a parent, regardless of the decay mode
The cell below selects Eu-152 and calls the dr_photon_total function which performs the task. Note that the function is called with a filter which cuts intensityes less than 2%
df = nl.pandas_df_nl(nl.nuclide("152EU").decays[0].dr_photon_tot("DR_PHOTON_TOTAL.INTENSITY > 2"), pd)
go.Figure().add_trace(go.Scatter(x=df['energy'], y = df['intensity'], mode = "markers")).update_layout(xaxis_title='Energy [keV]',
yaxis_title='Intensity [%]').show()
df
| parent_nucid | parent_l_seqno | energy | energy_unc | energy_limit | intensity | intensity_unc | intensity_limit | type | count | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 152EU | 0 | 6.3540 | 0.0000 | = | 14.000 | 0.500 | = | X | 1.0 |
| 1 | 152EU | 0 | 39.5220 | 0.0000 | = | 20.870 | 0.200 | = | X | 1.0 |
| 2 | 152EU | 0 | 40.1170 | 0.0000 | = | 37.800 | 0.300 | = | X | 1.0 |
| 3 | 152EU | 0 | 45.9980 | 0.0000 | = | 14.850 | 0.200 | = | X | 1.0 |
| 4 | 152EU | 0 | 121.7817 | 0.0003 | = | 28.530 | 0.160 | = | G | 1.0 |
| 5 | 152EU | 0 | 244.6974 | 0.0008 | = | 7.550 | 0.040 | = | G | 1.0 |
| 6 | 152EU | 0 | 344.2785 | 0.0012 | = | 26.590 | 0.200 | = | G | 1.0 |
| 7 | 152EU | 0 | 411.1165 | 0.0012 | = | 2.237 | 0.013 | = | G | 1.0 |
| 8 | 152EU | 0 | 443.9606 | 0.0016 | = | 2.827 | 0.014 | = | G | 1.0 |
| 9 | 152EU | 0 | 778.9045 | 0.0024 | = | 12.930 | 0.080 | = | G | 1.0 |
| 10 | 152EU | 0 | 867.3800 | 0.0300 | = | 4.230 | 0.030 | = | G | 1.0 |
| 11 | 152EU | 0 | 964.0570 | 0.0050 | = | 14.510 | 0.070 | = | G | 1.0 |
| 12 | 152EU | 0 | 1085.8370 | 0.0100 | = | 10.110 | 0.050 | = | G | 1.0 |
| 13 | 152EU | 0 | 1112.0760 | 0.0030 | = | 13.670 | 0.080 | = | G | 1.0 |
| 14 | 152EU | 0 | 1408.0130 | 0.0030 | = | 20.870 | 0.090 | = | G | 1.0 |
See that one can refer to a column in a dataframe by its position, in this way one has a template to plot 3d
# convenience function to plot 3D
def plot3d(fields, condition):
df = nl.pandas_df(fields,condition, pd)
fig = go.Figure(data=[
go.Scatter3d(
x=df[df.columns[0]], y = df[df.columns[1]], z = df[df.columns[2]],
mode='markers', marker=dict(size=8, color='blue') )
])
fig.update_layout(height=1000 ,scene=dict(zaxis=dict( type='log'))).show()
return df
fields = "LEVEL.NUC.Z , LEVEL.NUC.N , LEVEL.HALF_LIFE_SEC"
condition = "LEVEL.ENERGY = 0"
plot3d(fields, condition)
| z | n | half_life_sec | |
|---|---|---|---|
| 0 | 47 | 53 | 1.206000e+02 |
| 1 | 48 | 52 | 4.910000e+01 |
| 2 | 49 | 51 | 5.650000e+00 |
| 3 | 36 | 64 | 7.000000e-03 |
| 4 | 42 | 58 | 2.212141e+26 |
| ... | ... | ... | ... |
| 3349 | 92 | 122 | 5.200000e-04 |
| 3350 | 78 | 87 | 2.600000e-04 |
| 3351 | 80 | 90 | 8.000000e-05 |
| 3352 | 9 | 4 | 4.517205e-22 |
| 3353 | 17 | 35 | NaN |
3354 rows × 3 columns
fields = "DR_DELAYED.PARENT.Z , DR_DELAYED.PARENT.N , DR_DELAYED.PARENT_LEVEL.HALF_LIFE_SEC"
condition = "DR_DELAYED.TYPE = 'DN' and DR_DELAYED.PARENT.Z > 10 and DR_DELAYED.PARENT.Z < 100"
plot3d(fields, condition)
| z | n | half_life_sec | |
|---|---|---|---|
| 0 | 37 | 63 | 0.0520 |
| 1 | 19 | 34 | 0.0300 |
| 2 | 49 | 84 | 0.1650 |
| 3 | 51 | 84 | 1.6790 |
| 4 | 33 | 52 | 2.0210 |
| 5 | 19 | 30 | 1.2600 |
| 6 | 29 | 50 | 0.2410 |
| 7 | 31 | 52 | 0.3081 |
| 8 | 19 | 31 | 0.4720 |
| 9 | 12 | 21 | 0.0905 |
| 10 | 19 | 33 | 0.1100 |
| 11 | 19 | 32 | 0.3650 |
| 12 | 37 | 56 | 5.8400 |
| 13 | 37 | 59 | 0.2030 |
| 14 | 37 | 58 | 0.3777 |
| 15 | 37 | 57 | 2.7020 |
| 16 | 47 | 76 | 0.2990 |
| 17 | 53 | 86 | 2.2800 |
| 18 | 53 | 84 | 24.5000 |
| 19 | 53 | 85 | 6.2600 |
| 20 | 35 | 53 | 16.3400 |
| 21 | 17 | 29 | 0.2320 |
| 22 | 35 | 55 | 1.9100 |
| 23 | 11 | 19 | 0.0480 |
| 24 | 35 | 54 | 4.3570 |
| 25 | 35 | 52 | 55.6800 |
| 26 | 11 | 18 | 0.0441 |
| 27 | 11 | 21 | 0.0132 |
| 28 | 35 | 57 | 0.3140 |
| 29 | 37 | 60 | 0.1691 |
| 30 | 37 | 62 | 0.0540 |
For each nuclide, the properties daughters_chain and parents_chain contain the offsprings and the ancestors, respectively
The example shows the seamless interplay between sets of NDLab classes and pandas dataframes
# get the offsprings as NDLab class instances
dau = [n.gs for n in nl.nuclide("241AM").daughters_chain]
# dump in a dataframe
df = nl.pandas_df_nl(dau,pd)
# Plot. Hover on a point to see more info
px.scatter(df,
x= df.index,
y= df.half_life_sec,
hover_data=["nucid"],
labels={ "index": "Daughter # ", "y":"Half-life [s]"},
log_y=True).show()
fields = "LEVEL.J as j, LEVEL.ENERGY as e"
filter = "LEVEL.NUC.Z = 20 and LEVEL.JP_METHOD = RIPL_J_UNIQUE"
df = nl.pandas_df(fields, filter, pd).sort_values(by=['j'])
go.Figure().add_trace(go.Scatter(x=df['j'], y = df['e'], mode = "markers")).update_layout(xaxis_title='J',
yaxis_title='Energy').show()
fields = "GAMMA.START_LEVEL.J as j, GAMMA.ENERGY as e"
filter = " GAMMA.NUC_ID = '16O' "
df = nl.pandas_df(fields, filter, pd).sort_values(by=['j'])
go.Figure().add_trace(go.Scatter(x=df['j'], y = df['e'], mode = "markers")).update_layout(xaxis_title='J',
yaxis_title='Energy').show()
fields = "LEVEL.J as j, LEVEL.DIPOLE_MM as m"
filter = "LEVEL.NUC.Z % 2 = 1 and LEVEL.NUC.N % 2 = 0 and LEVEL.JP_METHOD = RIPL_J_UNIQUE"
df = nl.pandas_df(fields, filter, pd).sort_values(by=['j'])
go.Figure().add_trace(go.Scatter(x=df['j'], y = df['m'], mode = "markers")).update_layout(xaxis_title='J',
yaxis_title='Magn. dipole [mu N]').show()
fields = "LOG(L_DECAY.LEVEL.HALF_LIFE_SEC) as h , 1/ SQRT(L_DECAY.NUC.QA + L_DECAY.LEVEL.ENERGY) as e "
filter = "L_DECAY.NUC.Z = 92 and L_DECAY.NUC.N % 2 = 0 and L_DECAY.MODE = DECAY_A "
df = nl.pandas_df(fields, filter, pd).sort_values(by=['e'])
go.Figure().add_trace(go.Scatter(x=df['e'], y = df['h'], mode = "markers")).update_layout(xaxis_title='E(keV)<sup>-1/2</sup>',
yaxis_title='log(T<sub>1/2</sub>)').show()
# practice 1
DR_BETA.PARENT_LEVEL.JP
# practice 2
fields = " GAMMA.ENERGY "
filter = " GAMMA.START_LEVEL.JP = '2+' "
print(nl.csv_data(fields, filter))
print(nl.json_data(fields, filter))
# practice 3
fields = "GAMMA.NUC.Z as z, GAMMA.NUC.N as n , abs( GAMMA.MIXING_RATIO / GAMMA.ENERGY / 0.835 * 1000) as b"
filter = "( GAMMA.NUC.Z % 2 = 0 ) and ( GAMMA.NUC.N % 2 = 0) " # even-even nuclides
filter += " and GAMMA.START_LEVEL.JP = '2+' and GAMMA.START_LEVEL.JP_ORDER = 2 " # starts from Jp=2+ ,2nd occurrence
filter += " and GAMMA.END_LEVEL.JP = '2+' and GAMMA.END_LEVEL.JP_ORDER = 1 " # ends at Jp = 2+ ,1st occurrence
filter += " and (GAMMA.MULTIPOLARITY = 'E2+M1'or GAMMA.MULTIPOLARITY = 'M1+E2') " # E2 + M1 multipolarity
nl.csv_data(fields, filter)